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Abstract 

It is preliminarily known that Aedes mosquitoes are very close to 
humans and their dwellings, also give rises to a broad spectrum of 
diseases: dengue, yellow fever, chikungunya. In this paper, we explore 
a multi-age-class model for mosquito population secondarily classi¬ 
fied into indoor-outdoor dynamics. We accentuate a novel design for 
the model in which periodicity of the affecting time-varying environ¬ 
mental condition is taken into account. Application of the optimal 
control with collocated measure as apposed to the widely-used pro- 
totypic smooth time-continuous measure is also considered. Using 
two approaches: least-square and maximum likelihood, we estimate 
several involving undetermined parameters. We analyze the model 
enforceability to biological point of view such as existence, unique¬ 
ness, positivity and boundedness of solution trajectory, also existence 
and stability of (non)trivial periodic solution(s) by means of the basic 
mosquito offspring number. Some numerical tests are brought along 
at the rest of the paper as a compact realistic visualization of the 
model. 
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likelihood; optimal control; the basic mosquito offspring number. 


Mathematical Subject Classifications: 49J15; 49M05; 92D30; 

37B55. 

1 Introduction 

We consider a mathematical model of mosquito population dynamics within 
the framework similar to [28, 29], where the continuous population evolution 
is formulated as initial value problem (IVP) of non-autonomous system 

x = f(t,x, 11577), te[0,T],x(0)=x o 8 0. (1) 

In natural setting, the system equation contains a hyperparameter r/ whose 
appearance describes a collection of measurable intrinsic factors: e.g. natural 
births, natural deaths, age-based transitions. Studying the qualitative be¬ 
havior of the solution, one demands the fluctuation phenomena within those 
factors to be distinguished as those of which essentially depend on time. As 
a consequence of uncertainty in environmental condition, some elements of 
rj differ in time with possible seasonal trends: monotonic increasing, mono¬ 
tonic decreasing, oscillating, or even fluctuating with Brownian-type move¬ 
ment. Many references have even hypothesized that such intrinsic factors 
may behave with periodic streamline in many cases due to environmental 
changes: cf. [5, 11, 10, 12, 22, 23]. Ironically, in a national integrated 
mosquito management programme, for example, Bonds [3] summarized that 
fluctuating meteorology (raindrop, wind speed, air temperature, air humid¬ 
ity, terrestrial radiation, etc) had been out of concern during deployment of 
control devices in the field. Therefore, this costed substantial inefficiencies in 
mass disposal. This dependency of parameters in time brings the model into 
non-autonomous groundwork. As an extrinsic factor, the control measure u 
is incorporated into the system for which it plays as a system regulator to¬ 
wards achievement of the general objectives: i.e. minimizing both population 
and cost for the control. In other words, the following objective functional 



attains its minimum for given weighting constants {ca x ,i}ie/ x an d 
J x = {1, • • • , 5} and I u = {1, 2}. 

In line with matching the underlying dynamical process of the solution 
with that from empirical measurement, we demonstrate estimation of some 
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undetermined parameters in the model by firstly settling them as random 
variables. Then by exploiting information from the system, one can charac¬ 
terize the solution as a handling function expressed in terms of such random 
variables. We utilize the property of Maximum Likelihood (ML) and Least 
Square (LS) which is basically minimization of the difference between the 
handling function value-points and the measured data over all possible values 
of undetermined parameters in a bounded set. The next problem arises when 
analytical solution of the system, or eventually the handling function, can not 
be determined explicitly due to complexity of the equation. For this reason, 
a schematic Local Linearization (LL) method provides a trade-off between 
numerical accuracy and computational outlay. Recently, several regimes in 
coping with ML for parameter estimation within epidemic dynamical system 
have been explored e.g. in [31, 9, 6]. 

There have been numerous mathematical papers discussing application 
of optimal control scenario to the mosquito reduction issue (see e.g. [28, 
29, 8, 24, 7, 18, 30] and some references therein). The authors used proto- 
typic autonomous model utmost, encouraging us to propose a novel approach 
adopting non-autonomous dynamical system theory. In this paper we restrict 
our main scope to the application of temephos and ULV aerosol. Enhance¬ 
ment of indoor-outdoor dynamics and utilization of polynomial collocation 
design to the control measure are parts of our interests. The choice of the 
design aims at achieving minimization of the costly objective meanwhile pro¬ 
nouncing more efficient and accurate control deployment. 


2 Model and analysis 


In favor of IVP (1), x denotes the time-variant state which folds five con¬ 
secutive elements, each of which represents the number of: indoor eggs x 1; 
outdoor eggs x 2 , indoor larvae x 3 , outdoor larvae x 4 and adults x 5 . A control 
measure u is injected into the system as an active feedback regulator whose 
elements represent the impact rate for investment of: temephos u 4 and ULV 
aerosol u 2 . Note that we omit writing argument t when it is obvious. To go 
into details, the governing system (1) is unfolded as 


xi = pa(t)x5-ftxi-guiXi-//iXi, 

x 2 = (1 - p)a(t)x 5 - /3 2 x 2 - p 2 x 2 , 

x 3 = /Axi - 71 X 3 - /? 3 x 3 - uix 3 - ru 2 x 3 - /i 3 x 3 , 

x 4 = /3 2 x 2 - y 2 x 4 - /3 4 x 4 - su 2 x 4 - /i 4 x 4 , 

x 5 = /3 3 x 3 + /3 4 x 4 - u 2 x 5 - /i 5 x 5 . 


(3a) 

(3b) 

(3c) 

(3d) 

(3e) 
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All the involved parameters are positive and are briefly explained as fol¬ 
lows: p, q, r , s are the plausible probabilistic constants; a(t) = e + eo cos (at) 
(where e > e 0 > 0) is the birth rate of potential eggs depending qualitatively 
on meteorology distribution; /3{ 1 , 2 , 3 , 4 } are the age-transitional rates for the 
corresponding classes; 7 ( 1 , 2 } (where 71 > 72 ) are the driving forces to the 
arising competition amongst larvae and; ^{ 1 , 2 , 3 , 4 , 5 } are the death rates for the 
corresponding classes. 

It is assumed that the control measure u be in U '= C°([0, T]; B) a set of 
piecewise-continuous functions where B is a bounded block in R^. Since f is 
uniformly locally Liptschitz continuous on time-state domains and piecewise 
continuous on [0,T], then x should lie in C 1 ([0, T]; R 5 ) a set of piecewise- 

differentiable functions. It is defined a form £(f) = ||/||(i) by taking a 
norm ||-|| of / for each t. In this non-autonomous model, we initially denote 
x(£) = u(t, 0,x 0 ,u) the solution of (1) as a process. In a specified case 
when eo = 0 , typical analyses: existence and uniqueness, positive invariance, 
existence and stability of equilibria as regards the basic mosquito offspring 
number can be referred from our preceding work [29]. 

Consider A4 as a non-autonomous set where M. C [0 ,T] x R 5 . Denote by 
M t = f {x G R 5 : (t, x) G M} f-fiber of M for each t G [0, T], 

Definition 1 (Positively invariant) An autonomous set A4 C [0, T] x 

R 5 is called positively invariant under the process v if for any x e Ado, 
v(t, 0, x, u) G AA. t for all t > 0. 

Theorem 1 It holds that A4 = [0, T] x R'^ is positively invariant under the 
process f>. 

Proof Let n be (5 x 5)-matrix representing a collection of all normal vectors 
(by rows) to the boundary of R^_, dR+. It follows that n = —id where id is 
the identity matrix. Notice that at zth boundary, <9jR+, 

[nf (t, x, u)]^^ |UgI , < 0 uniformly Vi > 0 , i = 1 , ■ ■ • , 5. 

This generates evidence showing evolution of the solution-points in 9 r^_ in 
counter-direction or at least perpendicular to the corresponding normal vec¬ 
tor. Thus such trajectory of points, emanated from all t > 0, can not leave 
R®.. □ 

We next consider a decomposition over f as f (i, x, u) = A(t)x+cix|+c 2 x| 
to exemplify further analysis. In this case, A(t) is the Jacobian of the system 
evaluated at 0. Let V : [0, + 00 ) x C , 1 ([0, T]\ R^_) x C , 1 ([0, T\; R^_) —> R be a 
function defined as 

P(f,x, y)^(x-y) 2 . 
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This function delineates instantly three following conditions: (i) V > 0 if 
x/y and = 0 if x = y, (ii) V is uniformly locally Lipschitz continuous on 
dom(P) and (iii) for any {x n ,y n }„ 6N £ C' 1 ([0, T]; R®_), lim^,*, Vft, x n , y n ) = 
0 implies lim, woo £(x n - yj = 0 . 

Lemma 1 The following assertion applies 

lirn h~ 1 (V(t + h,x+hf( y t 1 x 1 u),y+hf(t,y,u))-V(t,x,y)) < g{t,V(t,x,y)) 

h-> 0- 

where g(t,w) is a continuous function exceeding 2 £{A(t))w for all t,w £ R+. 
In contrast, £(Aft)) denotes a corresponding form applied to a matrix Aft ) 
induced by if). 

Proof The inequality follows directly from straight-forward computation 
(ref. Cauchy-Schwarz inequality). □ 

Lemma 2 The following scalar non-autonomous equation 

w = gft,w) 

holds these two conditions: (i) gft, 0) = 0 for all t £ R + and (ii) for each r £ 
(0, Too), w = 0 is the only solution on [0, r] satisfying u>(0) = 0 beforehand. 

Remark In our model, A is a matrix-valued function over t whose elements 
contain the continuous function a and the piece-wise continuous function u. 
Along with the definition, it is clear that both a and u are bounded by some 
continuous function, i.e. there exists a continuous function / such that e.g. 
f?(u) < / for all t. Thus without lost of generality, 1 -form £\{A) induced by 
the 1 -norm ||-||i constitutes an easiness in proving that any induced form t{A) 
is bounded. This can immediately be seen by envisaging the matrix norm 
equivalence, in sense that there exists a bounded positive function Cft) such 
that £(A(t)) < C(t)£i(A(t)) for all t and any induced form £(■). Now the 
existence of such function g in the form g ft,w) = S(t)w can be drawn upon 
the fact that £(A) is bounded. 

We further state that Lemma 1 and Lemma 2 are inherently mild sub¬ 
stances to prove the uniqueness of solution in (1). The notion of uniqueness 
in non-autonomous dynamical system dates back to the seminal works by 
Murakami [15], Ricciardi-Tubaro [20] and Kato [14], where aforesaid condi¬ 
tions stipulated in Lemma 1 and Lemma 2 set the conforming requirements. 
Adopting materials in [13, 14, 27], we summarize the existence and unique¬ 
ness of the solution through the following corollary. 
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Corollary 1 (Uniqueness) IVP ( 1 ) has a unique solution x defined on R + 
which maps to R®_ for every Xq e R+. 

It is further imposed the following assumptions on (1): 

(HI) the hyperparameter r) is chosen to lie in the following set 

{<1 e R+ : x s (f;>j) > x 5 («;i;), x 4 (i;rj) > x 5 («;i;)V( € R + and ||i;|| < 00 } 

(H2) there exist positive continuous functional $ and sufficiently large con¬ 
stant L such that p lies in 

{h e R+ : 2 ^^-^) ^ x( ^ ; ^)) 2 - L + e R+j 


where c(t) = rna x{pa(t) — — 2 /q, (1 — p)ot{t) — @2 — 2 /i 2 , /h — — 

2/7-3, @2 Pa 2 /i 4 , a(£) + ^3 + /3 4 — 2 /r 5 }. 

Theorem 2 (Lyapunov method) Consider IVP ( 1 ). If there exist a Lya¬ 
punov function 'h : R + x C 1 ([0,T]; R®) — > R +; positive integers <^{ 1 , 2 , 3 }? posi¬ 
tive and nonnegative real numbers L, ( and positive continuous t-functionals 
${ 1 , 2 , 3 } where $1 is nondecreasing satisfying the following conditions 

(Al) di^xfi 1 < ^t(t,x) < d-fifixy 1 

(A 2 ) 4t(t, x) < —d 3 I(xY 3 + L 

(A 3 ) ^(t, x) — \l/(f, xfi 3 ^ 2 < ( 

for all f 6 R + , then the solution x is bounded on R + . 


Proof We aim at computing the following total derivative T\I /(t, x) exp(at) 
where a is a positive constant needed to be determined later. It is clear 

that ^\P(t, x) exp(at) = (\&(t, x) + a\l/(f, x)) exp(af) < (—$ 3 ifxfi 3 + L+ 
(Al) 

a'h(/, x)) exp (at) < 


^3 


a() exp (at) by taking a = inf fgR+ This results in 


. \ (A3) 

^(t, x) Q / 1 ' 2 + L + a\P(t, x) j exp (at) < (L + 


\l/(f, x) = ( \l/( 0 , x 0 ) + -— ^ exp (at) — -— + ^ ) exp (—at) 


(Ai) (L + aC) 

< $ 2 ( 0 )||x 0 |r 2 exp(-at) + 1 + U < $ 2 ( 0 )||xo| 


(L + a() 
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Using the left-side inequality in (Al), we get 


*(x) < ^i(O)- 1 ^ ^ 2 ( 0 )||xo||« + (L ^ aC) ^) 

since II \ is nondecreasing. □ 

Corollary 2 If (HI) and (H 2) are satisfied, then all solutions of (1) are 
bounded on R + . 

Proof This is a direct consequence of applying Theorem 2 to (1) with the 
following input: <qi, 2 , 3 } = 2, 'djig} = 1 , $3 = d, ( = 0 and \P(£, x) = £(x) 2 . □ 

Remark Now we have an appropriate bound for the solution namely ( ||xq|| 2 + L ) 2 . 


Definition 2 (Fundamental matrix) Consider a non-autonomous linear 
system 

z = W(t)z, t G R + , z(0) = z 0 (4) 

with Z(t, 0) as fundamental matrix for W(t), i.e. z(t,0,zo) = Z(t,0)zo, 
satisfying Z( 0,0) = id. It can further be verified that: (i) Z(t,s)Z(s, 0) = 
Z(t,Q), (ii) Z(t, O)^ 1 = Z(0,t) and (Hi) Z(t,t ) = id. 

Theorem 3 Consider (4) where W(t) is 2ir/a -periodic. There exist a dif¬ 
ferentiable 271/cr-periodic matrix Ti(t) and a constant matrix ^ such that the 
according fundamental matrix 

Z(t, 0) = Ti(t)exp(r 2 £). 

Remark Several highlighting insights regarding Theorem 3: 

• Ti(t),r 2 need not to be unique and real even though W(t) is real 

• the theorem holds for W ( t ) complex 

• since Z(0, 0) = id, or eventually rx(0) = T 1 (27r/cr) = id, then Z(2ti/ 
cr, 0) = Ti(27r/(j) exp(r 2 27r/cr) = exp(r 2 27r/cr). 

Definition 3 (Floquet exponent and Floquet multiplier) Let p(A ) be 

the spectrum of A. In Theorem 3, each element of p(T 2 ) is called as Floquet 
exponent meanwhile each element of p(Z(2ir/a, 0)) is called as Floquet mul¬ 
tiplier. 
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Take a look back to (1). Irrespective to the appearance of the control, 
f(f, x, 0) can be identified and decomposed as f(f, x, e 0 ) = fo( x ) + e 0 fi(£, x) 
where f 0 is autonomous and f i (t, x) = f* ^f(f, x, £e 0 ) d£. It is preliminarily 
known the following: (i) f (t + 27t/cx, x, eo) = f(t, x, eo) for all t E R+ and (ii) 
if Q is equilibrium point of x = fo( x ) then f 0 (Q) = 0 . 

Theorem 4 (Existence of periodic solution) Let f(t,x,e 0 ) = / 0 («c) + 

eo/i(t, a:) and Q be an equilibrium point ofx = f 0 (x). If 27i/ap(y x f 0 (Q)) E5 
2n/cr\ 2ixiZ, then there exist a neighborhood U(Q) and e± such that for ev¬ 
ery |e 0 | < ei, there exists a 2 tt/ a-periodic solution of ( 1 ) with unique initial 
value Xq = a%(e 0 ) EU(Q). 

Proof Let x(i, x 0 ,e 0 ) be such solution. Letting x 0 E Wi(Q) C R+ then x 
exists and is unique for all t E R + by Corollary 1 without lost of generality 

for |e 0 | < e^i. Let ((t) = V X 0 x(i, Q,0). Since V X 0 x(f, x 0 , e 0 ) deductively 
satishes the variational equation 

-j“V Xo x(f, x 0 , e 0 ) = V XQ f (£, x(t, x 0 , e 0 ), e 0 ) 

= V x f(t, x(t, x 0 , e 0 ), e 0 )V Xo x(t, x 0 , e 0 ), 

with V Xo x(0, x 0 , eo) = id, then ( immediately satisfies 

4c = V x f 0 (Q)C with C( 0 ) = id. 

at 

Let ^(xq^o) = x(27r/o-,x 0 ,e 0 )-x 0 . Since f E C\r + x R^_ x (-ei,i, eyi); R 5 ), 
then S E C 1 (R^_x(—ei,i, ei,i); R 5 ). It is clear that S( Q, 0) = 0 and V Xo S'(Q, 0) 
exp(V x fo(Q)27r/cr) — id. Let v be eigenvector of V x fo(Q) associated with 
eigenvalue A. We know that (exp(V x f 0 (Q)27r/cr) — id)v = (exp(27r/aA) — l)v 
making a clearance that exp(27r/cxA) — 1 is eigenvalue of V X 0 S'(Q, 0). If 2n/ 
crA ^ 27riZ then det(V X 0 S'(Q, 0)) = =1 exp(27r/crAj) -1^0, making 

V Xo S'(Q, 0) invertible. By implicit function theorem, there exist a domain 
W 2 (Q) x (— 61 , 2 , 61 , 2 ) and a smooth x 0 (eo) for (eo,x 0 (eo)) defined on this do¬ 
main such that S(xo(eo),eo) = 0 or eventually x(27t/ct, x 0 (eo), eo) = x 0 (eo). 
Since f is 27r/<r-periodic over t, then x(t + 27r/cx, x 0 (e 0 ), e 0 ) = x(£, x 0 (e 0 ), eo) if 

and only if x(2n/a, x 0 (eo), eo) = x 0 (eo). Now letting U(Q) == L/i(Q) r\U 2 (Q) 
and 61 such that (— 61 , 61 .) C (— 61 , 1 , 61 , 1 ) D (— 61 , 2 , 61 , 2 ) follows the desired 
domain of existence. □ 

Let bi = ep, b 2 = e(l - p), b 3 = /3 1 , b 4 = (3 2 , b 5 = (3 3 , b 6 = j3 A and 
di = (3i+ pi, d 2 = (3 2 + /i 2 , d 3 = f3 3 + fi 3 , d 4 = (3 4 + /x 4 , d 5 = p 5 . With the 


same technical arrangement using the next generation method [25] as in [29], 
we define 


n(d; h d 4 ) = 


b 1 hh + b 1 b ± h 

d§ d\ d% d§ d% d 4 


(5) 


as the so-called basic mosquito offspring number. In the domain of interest 
R+, it has been proved in [29] that two equilibria of x = f 0 (x) exist: zero 
equilibrium and a positive equilibrium Q. 


Lemma 3 There exist two 2 tt/ a-periodic solutions of (1) in R + : trivial so¬ 
lution x = 0 if 1 Z(d 3 , d.f) f 1 and nontrivial solution x = v associated with 
nontrivial autonomous equilibrium Q ifdZ(d 3 + 27 !* 3 , d 4 + 272 ^ 4 ) f 1 where 
f,-,xfxf-) = Q. 


Definition 4 (Stable periodic solution) A periodic solution v is said to 
be stable if for every e > 0 there exists S > 0 such that ||a^ — i/(r)|| < 5 implies 
I(x(t,r, xf) — u(t)) < e for all t >r. Additionally, if r, xf) — 

v{t)) = 0 then v is said to be asymptotically stable. 


Theorem 5 (Stability of periodic solution) Let y = x — v and irre- 

def def 

spective to the control, h(t, y ) = f(t, y+v)— f(t , u) — W(t)y where W(t ) = 
Va;/(f, u). If the following conditions hold: 

(Bl) 1(h) < KI(y) 2 for some constant K 

(B2) all Floquet exponents y of u correspond to W ( t ) lie off the open left- 
half plane in C (or the Floquet multipliers lie off the open unit disk in 
C) 

(B3) in the decomposition exp(J Q < W (s) ds) = Ti(£) exp(r 2 t), T 1 is bounded 
w.r.t. If) 

then v is asymptotically stable. 


Proof It follows that y satisfies 

y = W(t)y + h(t,y), 

where h (t, 0 ) = 0 and V y h(£, 0 ) = 0 . 

Observe that I(Z(t, 0)) = I fexp(J Q W(s) ds)J < £(Ti(t))I (exp(r 2 £)) < 

(B3) 

£(r 1 (£))O 0 exp(— \t) < Cexp(-Xt) for some positive constants C and A. 
Fix 5 where CK8 — A < 0. Define 0 < a == A — CK8 and b = f If 
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||y 0 1 | < b then by continuity of the vector field there exists r such that y ex¬ 
ists on [0, r] satisfying £{y{t)) < 5 for all t 6 [0, r]. Generating the solution, 

we obtain £{y) = £ \Z(t, 0)y 0 + Z(t, 0) f Q Z(s, 0) _ 1 h(s, y) dsj < [||y 0 ||C + 

CK5 f* exp(\s)£(y) ds] exp(— Xt). Employing Gronwall’s Lemma, we get 
exp(Af)£(y) < ||y 0 ||C , exp(C , A^t) or £{y) < ||y 0 ||C'exp(-a£) < <5 exp (-at) 
such that Hindoo f?(y) = 0. Using some extension theorem, it can be proved 
that there exists e such that y is defined on [ 0 , r + e) by continuity of £{y). 

For all points, without lost of generality, {t n }nGN\oo = {t+ (1 — l/n)e} ngN \oo, 
note that ||(y(U)|| < <5exp(— atf) < 6 for ||y 0 || < b. This contradicts maxi¬ 
mally of r. □ 


Lemma 4 The trivial periodic solution v = 0 of (1) is asymptotically stable 
whenever IZ{d 3 , d^) < 1. 


Proof Using direct computation, it can be shown that 


K = 

0 







(6) 



r 1 0 

0 0 

P £0 

sin (a 

t) 1 






0 1 

0 0 

(1 - p)eo sin 

{at) 




Ti (t) = 


0 0 

1 0 


0 



bounded w.r.t. £{■) 

(7) 



0 0 

0 1 


0 







. 0 0 

0 0 


1 







—d\ 

0 

0 

0 

bi 






0 

—d 2 

0 

0 

^2 




r 2 = 


63 

0 

—d 3 

0 

0 



(8) 



0 

&4 

0 

—c?4 

0 






0 

0 

bs 

h 

—d 3 





resulting r 1 (27r/cr) = id. Noticing [29] Theorem 3.2, it has been proved 
whenever IZ{d 3) d A ) < 1 then p(T 2 ) lies off the open left-half plane in C. □ 


Theorem 6 Let v be nontrivial periodic solution, ^{ 3 , 4 }(t) = \ ^{ 3 , 4 } (s) ds 

def 

where u {3A} (0) = ^{ 30 , 40 } and = mm te[ 0 i 2 ff /«T] ^{ 3 , 4 } (i). If 
IZ(d 3 + 2 7 l z/ 3 min , d 4 + 2 72 < n ) < 1 
then v is asymptotically stable. 
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Proof Let y = x — v. It is clear that y follows y = g(t,y) where 
g(t,y) = A(t)y + c 1 (yl + 2u 3 y 3 ) + c 2 (yl + 2u 4 y i ) by recalling our decomposi¬ 
tion upon f in (1). Now we can state that y delineates non-autonomous 
system with perturbance ^{ 2 , 3 } where 0 is trivial periodic solution. We 
can obtain easily the correspondence matrix for linearized system W(t) = 
V x g(i, 0 ). To save space, we briefly state that the fundamental matrix Z(t) = 
Z(t, 0) = exp(J Q W(s) ds) can not be presented easily as A^i) exp(A 2 i). 

To continue proceeding, the idea is by choosing Ai(£) c = Z(t) exp(—A 2 £) 
and A 2 for which exp(A 2 27r/<7) = Z (O) -1 Z (2n / a). Then this choice sat¬ 
isfies Z(t)Z(0)- 1 Z(2 tt/ a) = Ax(t) exp(A 2 t) • exp(A 2 27r/(T) = Ai(t + 2n/ 
cr) exp(A 2 (f + 2tt/ a)) = Z(t + 2n/a) which is nothing but the so-called Flo- 
quet theorem. Consider T 2 = ^(^ 3 ,^ 4 ) as in ( 6 ), it is clear that Z(27r/ 
a) = exp (2tt / aV 2 (d 3 + 2y\v 3 (2tt/ cr), d 4 + 2y 2 v 4 {2n/a))). We immediately 
get A 2 = T 2 (d 3 + 271 ^ 3 ( 2 ^/ 0 -), d 4 + 2y 2 v 4 (2'K/a)). If 7 Z(d 3 + 2y\v™ m ,d 4 + 
2 72 ^ 4 min ) < 1 then 7 ^(d 3 + 27 1 i/ 3 ( 27 r/cr), d 4 + 27 2 i 24 ( 27 r/cr)) < 1 , then with sim¬ 
ilar consequence as in Lemma 4, p(A 2 ) lies off the open left-half plane in C. 
Simultaneously, since 7 Z(d 3 + 2yiv 3 (t), d 4 + 2y 2 u 4 (t)) < 1Z(d 3 + 274 //™“, d 4 + 
2 l2 AD < 1 then Z is bounded and therefore Ai is bounded. The choice of 
zTW in [0, 27r/cr] by the fact that ^{ 3 , 4 } have the greatest deviations on their 
amplitude at this range. □ 

3 Parameter estimation 

Let rj be the hyperparameter of the model (1) and V C R 18 be its feasible 
region. Let X := G {1, • • • , 18} : p* unfixed} and 6 denotes a vector which 
collects all associated parameters whose indices in X. Let 0 C V respectively 
be feasible region for 6. The next key enabling technical simplification is that 
one can further rearrange the elements of rj as 77 = (r]J,6 T ) T . I 11 order to 
find an estimate of 9, it is essential to identify whether the system in nature 
is under control intervention or not. For the sake of simplicity, let us assume 
that there are no control treatments during matching process. Fixing r /f and 
setting u = 0 , we recast IVP (1) as 

±(t) = f(t,x(£);0), t G [0,T],x(0) = x 0 ,0 G 0. (9) 

Let J = f {0, • • • , N} and 

G n = {tj : tj = jAt, t N = tf, j G J} (10) 

be our set of discrete time-points. Taking a good solver for ODE, we assume 
that Eq. (9) results in the discrete process 0 :GatxOxR^_x0—*-r}_ mapping 
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(tj, 0, x 0 , 9) to Xj where the sequence conforms the regressing path. 

Let H : R+ —>■ R m be a function such that <f> = f Hocj) : x 0 x R+ x 0 —>■ R m . 

Assume (0,x 0 ) is fixed, leading to exposition $ : Gjy x 0 — > R m . 

Let /C C {1, • • • ,5}. Suppose that it is given a data set {T, x} of time- 
state points which folds the sample {T^, X)})^... k and let 

A = {i e 1C : 0 < t{ < T, l = 1 , • • ■ , hi}. (11) 

The first highlighting processes are given briefly as follows. In practice, since 
most T\ is beyond Gat, notorious interpolation and extrapolation processes 
are needed for all i e A and only interpolation process for all i e JC\A. The 
processes seek all corresponding state-points at all tj based on information 
from the known points given in the data set. Since an extrapolation process 
suffers from greater uncertainty, thus the higher ki for all i £ A will help to 
produce more meaningful results. Once we have the data set {T,x} ex- and 
interpolated with respect to aforesaid procedure, one draws the refined data 
set, {T, x}, where it holds {T* : j G J} = Gat for all i 6 1C and therefore 

X : Gat ->• rJ 1 . 

In contrast with easiness in real implementation and the lack of details 
in the data, we will always need such function H- like, H , which maps X 
from Ri|_ into R m as a collection course. We assume that the data for all 
state-classes in the model are not necessarily known. Now the corresponding 
number m should be taken to satisfy 1 < m < min{|/C|,5}. Working with 

def — 

the same treatment as in the regressing path, we let X = H o X : Gn —> R m , 
making X and $ comparable. 

3.1 Least-square approach 

Let us define the error of measurement 

jeX (12) 

Let e(9) = [e±(9) ■ ■ ■ e\j\(9)}. Given an estimate for 0, now our problem reads 
as 

find 9 e 0 such that J{9) = f ||e(6 ) )||p —> min. (13) 

In this formulation, ||-||i? denotes the Frobenius norm. 

3.2 Maximum likelihood approach 

Assume that {^j(9)}jej are considerably independent and identically dis¬ 
tributed (iid) since it most probably appears commonly that the data is 
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randomly distributed relative to the regressing path. Then we can assume 

{eMheJ ~ A/"(0, £) with the corresponding probability density function 
(pdf) 

<p( e 3> = fo ^ \ eX P ( 14 ) 

(27t)t det(S )2 L 2 J 

where A j = a/ ej(9) T Y~ 1 ej(9 ) is Mahalanobis distance from 6j(0) to 0 and 
det(£) is determinant of £. The joint pdf (jpdf) for all random variables 
is g iven b Y 


^M) = n (0 

(27r) 2 det(S )2 


exp 



1 

= , . m \J\ , x 1 J\ exp 

-hiAiii 

Z J 

(27r) 2 det(£) 2 

z 


(15) 


where A = (A 1; • • • , Au|) T . Independent from £, we get the fact 

||e(0)||ir —» 0 if and only if ||A|| 2 — > 0 if and only if ip —y max. (16) 

So the most sensible way of finding a good 9 is by maximizing <p, or by 
maximizing log<^, since log is monotonically increasing. 

Theorem 7 (Maximum likelihood for multivariate normal distribution) 

Let {tj(0)}j£j ~A/"(0, £) and 

def ' 

j&J 

then £(0) = argmax s log<^(e; 9) where <p(e]9) is given as in (15). 


so?) =' Y. m = Ts(») 


Remark The (m x m)-random variable S (9) is called as Wishart ma¬ 
trix which follows Wishart distribution W m (|j7"|,£) with parameters m (di¬ 
mension of the matrix), \J\ (degree of freedom) and £ (positively defined 
covariance matrix). We can always perform the standardization S (9) ~ 
W m (|J|,£) £-5S(0)£-^ ~ Wm(\J\,Im) Where W m (\J\,I m ) denotes the 

standard Wishart distribution. In special case when m — 1, Wi(|j7’|,l) = 

tfjv 

Hence the final constrained optimization problem reads as: find 9 G 0 
such that 


J(9) = f log </?(e; 6) = log 


(27r)^det(£(0))^ ^ 


|||A||! 


—>■ max. (17) 


In this case, A = (Ai, • • • , A^) 7 and A j = J £j(9) T Y,(9) l £j{9). 
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Remark In line with computation of optimal solution using a derivative- 
use method, one has to find the so-called Fisher’s score function f( 0) which 
is nothing but the Jacobian V@ log</?(e; 6) and (for Newton/quasi-Newton) 
the information matrix T\j\(0) which is negative of Hessian log</?(e;0). 
These computations require very lengthy expression and therefore one has 
to achieve very expensive evaluations. Nevertheless, heuristics should offer 
trade-off in direct solving but limit their speed in convergence. Another im¬ 
portant aspect in this problem is that, by giving 0 from the scratch, the 
value of the parameter 0 on each iterate seems converging to the boundary 
of 0. Initiatively, in this paper, we impose fixed value for all parameters right 
up in front, i.e. 77 = ( r/J,0j ) T , and then perturb the resulting solution with 
Gaussian noise along with the covariance matrix Ey. Matching the original 
with this perturbed model, one can perceive the process as Of- recovery. Con¬ 
sidering the estimate for 0 , there would be 2 possible methods which can be 
used: Wald confidence and the profile likelihood confidence methods. For \J\ 
very large, the variance var {Of) rsj Yj\(0f) where [var( 0 /)] ifc = [cov( 0/ ii5 0 /,*)], 
i,k = 1, ,|X|. Let Z\- T be (1 — r)-quantile of a standard normal dis¬ 

tribution and diagZ be a vector composed by selecting out main diagonal 
elements of Z. Choosing appropriate r, we gain Wald confidence interval 
running from Wald test: (H 0 : 0 = Of vs Hi : 0 7 ^ Of) as 


0 = 



Zi-Amgjfffoi), + 21 -rdiagy 



(18) 


One thing we need to make sure that at this large \ff \, numerical evaluation 
of the inverse information matrix should not be really expensive - one can 
approach it with numerical approximation on derivatives. Another method 
which is more accurate than Wald confidence method for \J\ small is the 
profile likelihood confidence method. The profile likelihood confidence in¬ 
terval (also called the likelihood ratio confidence interval) derives from the 
asymptotic Chi-square distribution of the likelihood ratio statistics. Let Z(@) 
and u(Q) denote the lower and upper bound of 0 respectively. It is known 


2 log 


( 9f) \ 

V ) 


< Am— 1;1—r 


which essentially determines 


(19) 


Z(0) = arg 


|(^(e; 0) : <p(e; 0) = </?(e; 0 f ) exp 


WV I—1;1 



( 20 ) 


For the upper bound u(Q), we take some arbitrary value in (Of, 00 ). To 
counteract the solution reaching the boundary of the given set, we estimate 
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a small positive number £ and refine the objective in both (13) and (17) using 
interior point function as: find 9 G 0 such that 

J{9) = ( ||e(6 l )||| - e [log(6 l - /(©)) + log(w(0) - 9 )] -»■ min (LS) 

and analogously: find 6 G 0 such that 

J(&) = f log <^(e; 9) + e [log(6 l — /(©)) + log(w(0) — 6 1 )] — > max. (MLE) 


3.3 Evaluation of $ 

In order to evaluate $, we adopt the property of Local Linearization (LL) 
method as it persuades balance between computational outlay and conver¬ 
gence issues. Related to as in [2, 17], the authors suggested to find the 
solution of 

x(i;0) = f(tj,x(tj-;0)) + V x f(^-,x(t,-;0))(x(t;0) -x(t,;0)) (21) 

on each subinterval [tj,tj+ 1 ) where x(£ 0 ) = x 0 and tj,tj+i G Gjv- The solution 
of (21) is given as the following recursion 

Xj+i(0) = Xj(0) + (exp [V x fj(0)At] - id) V x fj(0) -1 fj(0) (22) 

providing that V x fj(0) is invertible. In this formulation, x ? (0), f j{9) and 
V x fj(0) are abbreviations for x(t J ;6 l ), f(tj, x(ty 6 1 )) and V x f (tj, x(^-; 9)) re¬ 
spectively. 

Lemma 5 Let us devote to two solutions on the subinterval [tj,tj+ 1 ). Let 
4> be a process representing the analytic solution of (9) and (pAt be a pro¬ 
cess generated as the solution of (9) using LL method. Assume that f is 
uniformly Lipschitz continuous over all prescribed domains of its arguments. 
Independent from 9, there exists a positive constant C such that 

110 - ^AtlURto.tj+O) < CAt 2 uniformly Vtj,t j+1 G G N . 


Remark This Lemma exhibits the evidence that the smaller At taken in 
numerical computation, the more solution from LL method tends to analytic 
solution. One can take a look for the analogous proof of this Lemma in 
e.g. [19]. Another important problem needed to be tackled is how we can 
efficiently compute matrix exponential in (22). Interested reader can take a 
look into Pade approximate for matrix exponential, see e.g. [1, 26]. However, 
in this paper we omit the details of this approximate. 
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4 Optimal control problem 


4.1 Polynomial collocation 

Let us assume that the control measures are applied in every n days. The 
spacing time between times of application is assigned as h+ne where h, e G R 2 
and e is a vector containing unities. Let r 1,fc == /^(h + ne) and r 2,k = 
/fc(h + ne) + ne be two discrete time-points where Ik is a (2 x 2 )-diagonal 
matrix containing counter. Let * and */ denote respectively the MATLAB 
pointwise multiplication and division between two vectors. If diag(Jfc) counts 
all elements of the set {Oe, le, 2 e, • • • , (T — n) * /(h + ne)} in a consecutive 
manner, then both r 1,k and r 2,fc count some distinct numbers in R 2 . For the 
sake of simplicity, assume that (T — n) is divisible by h + ne in correspond 
to the operator */. Therefore, we have a finite collection of intervals 


def j- r gfc 2,k 


T-OH r r 

T = Ifr 


i 


Ti 


£ ) x [r. 


l,k 2 ,k 

To 


C )} 


fce{0,--- ,\\{T—n)*/ (h+ne)||oo} 


c Rj 


(23) 


Note that for all k G {min ((T — n) * /( h + ne )), • • • , ||(T—n)*/(h+ne)|| 00 }, 
it appears either [r, 1 ^, rf' k ) = 0 or [t^jT^') = 0 simultaneously since one 
may need h containing distinct elements. A collection of such counters should 
have zero cardinality. Let 8 ^ 2 (0 be some vector-valued function following 


def 


PM (MjGtrJVi’VthVf) 

0, (t,t) G ([0, T\ x [0, T])\ fe’Vi’*) x 


Ki-At) = 

y W, K u, U) ^ U w , - 1 j ^ L w , ^ \) \ ^L'l ) '1 ) ^ L'2 ) '2 J 

'(24) 

where p fe (£) is a vector-valued polynomial of degree being arbitrary. Let 
Sri-At) = 0 T ", Tl ((),--- ,CA* /(h+ " e) "“M) and 


P(t) = (Po(t)r-‘ >P||(t — n)*/(h+ne)||oo W) 

be (2 x [\\(T — n) */ (h + ne) ||oo +1])-vector-valued functions. Thus we design 
our control measure as 

c (t) '= v * 8 T i. T 2 (t). (25) 

In this case, v denotes a vector containing control measure values at T whose 
dimension is similar to that of p(i). Given a continuous-time control u(i). 
Dehne a p-collocation A(-; p) : U —> U such that for c = A(u; p), there is a 
weighting vector v satisfying c (t) = v * 5 T i. r 2(£) and ||u — c|| —y min for all 

{t, t) g r. 
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4.2 Existence of optimal control 

Designate the transformation over time-state variables on the following per¬ 
formance 

Xj i —y y { for all i = 1, • • • , 5 and t )->■ y 6 . 

As a consequence, there exists a function Y such that the non-autonomous 
equation (1) is similar to the following autonomous equation 

y = Y(y, u), yj(0) = x 0 for all i = 1, ■ ■ ■ ,5 and y 6 (0) = 0, t G [0, T] 

(26) 

where it holds Y 6 (y, u) = 1. Define V = {y : y = Y(y, u), t G [0, T], y(0) = 
y 0 , u G U} as a set of admissible states. Recall our objective functional 
J(u) = / 0 T g(y, u) dt as in (2) along with this transformation and compose 
the optimal control problem as 

find (y, u) G V x U such that J(u) — » min . (OC) 

The following lemma derives one appropriate material to prove existence of 
optimal control in (OC). 

Lemma 6 The following set 

S(t, y ) = {g(y , u) + 7 , Y x {y, «),■■■, Y 6 (y, u) : 7 > 0, u G B} 
is convex for all (t, y) G [0, T] x R® . 

Proof Fix (t, y) G [0, T} x as an arbitrary choice and write g(y, u) +7 = 
/c(ui, u 2 , 7 ). Keeping in mind that Y — Y(y, Ui, u 2 ) and n = k(ui, u 2 , 7 ) are 
continuous over u and 7 . Let B = [0, cq] x [0, a 2 ] such that (u l5 u 2 ) G B. Now 
consider that S(t, y) is a set of points (GR' where it structure can be stud¬ 
ied as follows. For fixed u 2 = 0 and 7 = 0, it is clear that £1 = k(ui,0,0), 
£{2,4} = K| 2 j 4 }(y, ui, 0) and £{3,5,67} are constant. This means that such 
points generate a parametric curve in R' whose projection on each £i£ 2 - 
and £i£ 4 -plane are convex quadratic, meanwhile on each £i£ 3 -, £!£ 5 -plane 
and so on are straight segments since [ 0 , af\ is bounded. If 7 goes from 0 
to 00 then this convex curve moves along £i-axis from initial position to 
infinity. At this stage, the generated 2 D-hyperplane, say P 2 d, is clearly con¬ 
vex. Moreover, for constant £{3,5,67} we can identify P 2D in £ 1 £ 2 £ 4 -Cartessian 
space. If u 2 goes from 0 to a 2 , then ? 2 D simultaneously moves along new 
axes called £5- and £ 6 -axis. It is clear that the set P 3 d = {(£1,£2,£4,£5) £ 
R 4 : (£1,£2,£4) e P 2D ,£5 e [Y 4 (y,Ui,a 2 ), Y 4 (y,Ui,0)]} is convex, and there¬ 
fore is the set P 4D = {(£1, £2, £4, £5, £e) e R 5 : ( 6 , 6 ,£5) e P 3 D,£e e 
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[y 5 (y,Ui,a2),y 5 (y,u 1 ,0)]}. Then for fixed (t, y), the set S(t, y) = {( 6 R' : 

( 6 , 6 ,£ 5 ,fe) e P 4 D,6,£ 7 constant} is also convex. □ 

Now we are ready to prove the existence of optimal control for our model. 

Lemma 7 There exists the only optimal pair (y, u) for the optimal control 
problem (OC). 

Proof We refer to Filippov-Cesari’s Theorem [4] to prove the existence of 
optimal pair, ft states whenever the following conditions hold 

1. There exists an admissible pair 

2. The set S(t, y) defined in Lemma 6 is convex for every (t, y) e [0, T] XR 6 

3. C/ is compact 

4. There exists a number 5 > 0 such that every solution £(y) < 8 for all 
t G [0,T] and all admissible pairs (y,u), 

then there exists such optimal pair. Ads 1 and 3 are trivial, meanwhile ad 2 
is proved in Lemma 6. Clearly a well-defined vector field Y in (26) conduces 
continuity of y on [0, T]. By Bounded Value Theorem, one can easily show 
that y is bounded on the compact f-domain [0, T], This completes the proof. 

□ 

Definition 5 (Saturation) Let (7([0, T]; R+) denote a set of piecewise-continuous 
functions mapping [0, T] into R+. It is defined the saturation 

T : C([0, T]; R^_) ->■ U = f C([0, T]; B) 

where the block B — [0, ai] x [0, a 2 ] as 

T(u) = max (0, min (a, u )). (27) 

Lemma 8 Let u * e C([0, T]; R^_) be the optimal control of broadening prob¬ 
lem (OC) by expanding the space for admissible controls. Let 

T:C([0,T]-R 2 + )^U 

be a saturation over the control resulting a projected optimal control u G U, 
or u = T (u*). Then (y, u) G T> x U is the solution of the original problem 
(OC). 

Working with the similar technical arrangement as in [29] Sec. 4, we generate 
an appropriate necessary condition for the optimal control as follows. 
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Theorem 8 (Necessary condition) Consider the broadening optimal con¬ 
trol problem (OC) by expanding the space for admissible controls. Let u * £ 
C([0, T];R+) be the minimizer for J and y * £ (7 1 ([0, T]; R®) be the resulting 
state. There exists a dual variable z k £ C <1 ([0, T]; R 6 ) such that the tuple 
(y*, z *, u*) satisfies the following system 


y = 


z = 


u = 


on 

dz 


with y*(0 ) = y 0 h 0, 


( y*,z*,u*) 

on 

dy 


( y*,z*,u *) 

( dU \ 

arg zero I —— I 


, z*(T) = 0 


(»*,•*) 


(28) 

(29) 

(30) 


def 

for all t £ [0, T], The functional TL(y, z, u) = y(y, u) + ( z , 7r(y, u)) zs Hamil¬ 
tonian functional, meanwhile all equations in (28) are respectively the state, 
adjoint, gradient equations and transversality condition. 

The adjoint (with transversality condition) and gradient equations can now 
be unfolded as 


Zi = 

-w x ,iyi + {/3i + qui + yi)zi - / 3 iz 3 , 


zi(T) = 1 

0 (31a) 

z 2 = 

— w x,2y 2 + (02 + t t 2)z 2 — / 3 2 Z 4 , 


z 2 (T) = 1 

0 (31b) 

Z 3 = 

-Wx, 3 y 3 + (27iy 3 + As + Ui + ru 2 + , 

q 3 )z 3 - /3 3 z 5 , 

Z3CO = 1 

0 (31c) 

Z -4 = 

— ^X,4y4 + (272y 4 + A 4 + SU 2 + /i4)z 4 

— At z 5, 

z 4 (T) = 1 

0 (31d) 

z 5 = 

—w x,5y5 + ( u 2 + ^ 5)25 




— 

[e + e 0 cos(cry 6 )]pzi - [e + e 0 cos(ay 6 )] (1 - p) z 2 , 

z 5 (T) = 1 

0 (31e) 

Z 6 = 

o-e 0 sin(cry 6 )py 1 Zi + < 7 e 0 sin(oy 6 )(l - 

p) y 2 z 2, 

z 6 (T) = - 

0 (31f) 

and 






w u ,iUi - gy 3 zi - 

y 3 z s = 0 


(32a) 


w u ,2U 2 - ry 3 z 3 - sy 4 z 4 - 

y 5 z s = 0 . 


(32b) 


The following algorithm illustrates our scheme to solve (OC). 


Remark Note that a termination criterion is included in the Line 11 namely: 
when A < Aq for sufficiently small Aq then the program terminates. 
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Algorithm 1 Gradient-based method for solving (OC). 

Require: k = 0, an initial guess for the control u fc G U, an error tolerance 
e > 0 and an initial step-length A > 0. 

1 : Compute the p-collocation u k <— A(u fc ;p). 

2: Compute yjj t— y fc (-;u fc ) and z k u t— z fc (-; y fc (-; u fe )) consecutively from 
the state (with forward scheme) and adjoint equation (with backward 
scheme). 

3: Compute the objective functional J( u k ). 

4: Compute u* fc (y(j, Zy u ) from the gradient equation and set u* fc G- 
A(u* fc , p). 

5: Update u fc+1 (A) u fc + Au* fc and u fc+1 •<— Y(u fc+1 ). Compute y(j +1 and 

Z y,u • 

6: Compute J(u fc+1 ) and set A J i— J(u fc+1 ) — J(u fc ). 

7: if |AJ| < e then 

8 : Set (y, u, J) t— (y„ +1 , u fc+1 , J(u fc+1 )) then stop. 

9: end if 

10 : while A J > 0 do 

11 : Update new A <— arg min s6 [ 0A ] 0(s) := J(u fc+1 (s)) where 0 is a 

quadratic function as representative of J with respect to the step- 
length s. Note that the solution exists since 0(0), 0'(O) and 0(A) can 
be computed directly. Compute new u fc+1 (A) *— u k + Au* fc and set 
u fc+1 X( u fc+1 ). Then compute new y^ +1 and Zy+h 
12 : Compute J(u fc+1 ) and set A J J(u fc+1 ) — J(u fc ). 

13: if |AJ| < e then 

14: Set (y, u, J) <— (y(j +1 , u fc+1 , J(u fc+1 )) then stop. 

15: end if 

16: end while 

17: Set k <— k + 1. 

18: Go to Step 4. 

19: return The tuple (y, u, j). 


5 Numerical tests 

Table 1 gives an estimate value of all parameters used in the model. In trial 
scheme, we aim at recovering 2 parameters: 6 1 = e 0 and 0 2 = p. We use 
H(x) = H(x) = Xj, T,f = 10 uniformly for all classes. We run ge¬ 

netic algorithm as a core program to solve both (LS) and (MLE) with the 
following computer specification: operating system OSX 10.9.4, processor 
2.6 GHz Intel Core i5, RAM 16 GB, programming language Python 2.7.6 
64bits with programming environment Spyder 2.2.4 and accuracy 4 digits. 
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The corresponding result, can be looked np in Table 2. The role of two impor¬ 
tant parameters in the model to magnitude of the basic mosquito offspring 
number is shown in Fig. 1. Concurrently, from Figs. 2, 3 and 4 we show 
the performance and attainment of optimal control to suppress the size of 
mosquito population. 


e e 0 p 

q r 

s fix 

h 

ft 

ft 

7i 

72 

hi 

p 2 

3 2 0.4 

0.04 0.05 

0.05 0.3 

0.2 

0.08 

0.05 

0.004 

0.0026 

0.02 

0.01 


h3 

h 4 

h5 

T 

h Tl <n x ,{l,2,3,4} <^x,5 Wj,{1,2} 

a 

x 0 

<T 

0.02 

0.01 

0.4 

151 

[9,19] 1 1 2 4xl0 3 

[1,1] 

[21,43,24,37,8 

‘2n 

|r/4| 


Table 1: Estimate for all parameters involved int the model (1). 


Problem 


Scheme 1 

Scheme 2 

Scheme 3 

Scheme 4 

Scheme 5 


IT 

100 

1000 

2000 

5000 

10000 

(LS) 

ET 

17 s 

125 s 

323 s 

901 s 

1987 s 

At = 0.1 

Ox 

1.7213 

1.7326 

1.8001 

1.8995 

1.9021 


o 2 

0.2001 

0.2567 

0.3210 

0.3334 

0.3789 

(MLE) 

ET 

191 s 

295 s 

692 s 

1443 s 

3403 s 

At = 0.1 

9i 

1.777 

1.8743 

1.9123 

1.9375 

1.9786 


o 2 

0.2932 

0.3031 

0.3635 

0.3800 

0.3994 

(LS) 

ET 

185 s 

1367 s 

3715 s 

10101s 

27456 s 

At = 0.01 

9i 

1.8514 

1.9522 

1.9812 

1.9991 

1.9999 


02 

0.3465 

0.3821 

0.4187 

0.4097 

0.4091 

(MLE) 

ET 

2921 s 

3229 s 

7942 s 

19534 s 

48047 s 

At = 0.01 

0 i 

1.9987 

1.9998 

1.9999 

2.0000 

2.0000 


02 

0.3985 

0.4023 

0.4003 

0.4000 

0.4000 


Table 2: The standard genetic algorithm (GA) performance in solving (LS) 
and (MLE). IT = number of iterates, ET = elapsed time, s = second. 


6 Concluding remarks 

We have exhibited the mosquito population dynamics model using some uni¬ 
fying theories bearing from non-autonomous dynamical system. Imposing 
relevant assumptions over all parameters in the model, we prove positivity, 
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Figure 1: The basic mosquito off¬ 
spring number lZ(d^,di) in (p, e)- Figure 2: Optimal control that re¬ 
plane. suits to Fig. 4. 




Figure 3: Uncontrolled dynamics. 


Figure 4: Controlled dynamics. 
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uniqueness and boundedness of the corresponding solution. With the Floquet 
theory, we prove that trivial periodic solution exists if the basic mosquito off¬ 
spring number TZ(d 3 ,d 4 ) ^ 1 and is asymptotically stable if 7?.(<73, d 4 ) < 1. 
In this paper, we can not derive a direct relationship between existence 
and stability of nontrivial periodic solution v correspond to nontrivial au¬ 
tonomous equilibrium Q and the basic mosquito offspring number. It states 
whenever max {7 Z(d 3 + 27^3, d 4 + 272X4), 7 Z(d 3 + 27 4 z/™ n , d 4 + 27 2 z/™ m )} < 
1 (see Lemma 3 and Theorem 6) then v exists and is asymptotically stable. 
From Table 1 and Fig. 1, all set of parameters satisfying lZ(d 3 , d 4 ) > 1 makes 
nontrivial periodic solution attracting the original solution, therefore it is 
graphically asymptotically stable. From Fig. 1, one has to reduce e (meaning 
that one has to ensure that meteorology does no longer support mosquito 
life) in order to reduce the basic mosquito offspring number significantly. A 
condition when lZ(d 3 , <73) < 1 guarantees that the mosquito population can 
return in insignificant number within finite time and completely die out for 
expanded time. 

Parameter estimation in our paper bears from necessity of suitable codes 
which have extreme reliability in real implementation. Three generic meth¬ 
ods: Local Linearization (LL), Pade approximate and Genetic Algorithm 
(GA) come into the play of estimating some parameters in the model. Sum¬ 
marizing the performance of our codes (ref. Table 2), one shows that the 
compact program executes in exponential time with the low convergence on 
average. Meanwhile, it highlights in the table that MLE scheme converges 
more rapid than LS scheme with respect to the number of iterates. Beyond 
naive implementation of GA, we examine that this low convergence results 
from expensive evaluation of the objective, even rigorous computation of LL 
solution on each iterate. Nevertheless, heuristic GA has gained a good deal 
in popularity of solving optimization problem with expensive objective. We 
strive to work on finding efficient codes of derivative-use method as a striking 
resemblance with GA. This in turn enables us to compare both methods in 
pursuit of efficient codes within the same framework, even with the higher 
size of both the system and undetermined parameters. 

We present some brief conclusions from application of optimal control. In 
this paper, our work is circumscribed by application of constant collocation. 
This can be more or less a key step toward development of efficient deploy¬ 
ment of the control. Further application of polynomial collocation of degree 
> 1 is needed to propose well-suited program, being readable throughout 
academia. By comparing Figs. 3 and 4, we conclude that optimal control 
can generally reduce all classes in the mosquito population. From Fig. 2, 
it is noticed that application of the ULV aerosol has to be enhanced rather 
than that of temephos. It is concluded that the pattern of optimal control 
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fluctuates with same tendency as that of the model dynamics. This means 
that the dynamics pattern of the model has a strong influence on the control. 
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